
%%% Weighting the Evidence: A Rank-Dependent Model of Outdoor Recreation, June 2024
%%% This function evaluates the log-likelihood function and its gradient for the RDU utility model (Prelec weighting function)

function[LnL, gradient]=likelihoodGrad_RDUT_Prelec(K,Q,G,Z,P,Y,C,I,IJ,ID,parameters)

delta=[parameters(3), parameters(4), parameters(5)]';
vartheta=[parameters(9),parameters(10),parameters(11),parameters(12),parameters(13),parameters(14), parameters(15)];
u_W=zeros(length(Q),K); 
Wk=zeros(length(Q),K); 
cumP=cumsum(P,2,'reverse'); 
cumP=round(cumP,7); 

wP=exp(-(-log(cumP)).^parameters(16));
dwP=-wP.*(-log(cumP)).^parameters(16).*log(-log(cumP));

dwP(isnan(dwP) | isinf(dwP))=0;
piP=horzcat(-diff(wP,1,2),wP(:,end));
derivativeGamma=horzcat(-diff(dwP,1,2),dwP(:,end));

for k=1:K
   Wk(:,k)=parameters(1).*Q(:,k)+parameters(2).*G(:,k)+Z*delta;
   u_W(:,k)=(1-exp(-parameters(6).*Wk(:,k)))/parameters(6); 
end
F=sum(u_W.*piP,2)+parameters(7).*C; 
EF=exp(F);

f_theta1=sum(((1-parameters(6)*u_W).*Q).*piP,2);
f_theta2=sum(((1-parameters(6)*u_W).*G).*piP,2);
f_theta3=Z(:,1).*sum((1-parameters(6)*u_W).*piP,2);
f_theta4=Z(:,2).*sum((1-parameters(6)*u_W).*piP,2);
f_theta5=Z(:,3).*sum((1-parameters(6)*u_W).*piP,2);
derivativeAlpha=(((1/parameters(6))*(Wk+1/parameters(6))).*(1-parameters(6)*u_W))-1/(parameters(6)^2);
f_theta6=sum(derivativeAlpha.*piP,2);
f_theta16=sum(derivativeGamma.*u_W,2);

Eopt_out=exp(parameters(8).*ones(length(C),1)+Y*vartheta'+parameters(7).*C);

ef=EF.*(1-IJ);
eopt_out=Eopt_out.*IJ;
Denom=ef+eopt_out;
InactiveCosts=C.*IJ;
sumtheta1=f_theta1.*EF.*(1-IJ);
sumtheta2=f_theta2.*EF.*(1-IJ);
sumtheta3=f_theta3.*EF.*(1-IJ);
sumtheta4=f_theta4.*EF.*(1-IJ);
sumtheta5=f_theta5.*EF.*(1-IJ);
sumtheta6=f_theta6.*EF.*(1-IJ);
sumtheta7=C.*EF.*(1-IJ);
sumtheta16=f_theta16.*EF.*(1-IJ);

sume=accumarray(ID, Denom); 
sumeInactiveCost=accumarray(ID, InactiveCosts);
sumActiveTrips=accumarray(ID, ef);
sumInactiveTrips=accumarray(ID,eopt_out);
Sumtheta1=accumarray(ID,sumtheta1);
Sumtheta2=accumarray(ID,sumtheta2);
Sumtheta3=accumarray(ID,sumtheta3);
Sumtheta4=accumarray(ID,sumtheta4);
Sumtheta5=accumarray(ID,sumtheta5);
Sumtheta6=accumarray(ID,sumtheta6);
Sumtheta7=accumarray(ID,sumtheta7);
Sumtheta16=accumarray(ID,sumtheta16);
SUME=kron(sume,[1 1 1]');
SUMEInactiveCost=kron(sumeInactiveCost,[1 1 1]');q
SUMEActiveTrips=kron(sumActiveTrips,[1 1 1]');
SUMEInactiveTrips=kron(sumInactiveTrips,[1 1 1]');
SUMEtheta1=kron(Sumtheta1,[1 1 1]');
SUMEtheta2=kron(Sumtheta2,[1 1 1]');
SUMEtheta3=kron(Sumtheta3,[1 1 1]');
SUMEtheta4=kron(Sumtheta4,[1 1 1]');
SUMEtheta5=kron(Sumtheta5,[1 1 1]');
SUMEtheta6=kron(Sumtheta6,[1 1 1]');
SUMEtheta7=kron(Sumtheta7,[1 1 1]');
SUMEtheta16=kron(Sumtheta16,[1 1 1]');

F=ef./SUME; 
FTerm=F(F~=0);
FTermI=I(F~=0);
S=eopt_out./SUME; 
STerm=S(S~=0);
STermI=I(S~=0);

LnL=-(sum(log(FTerm).*FTermI)+sum(log(STerm).*STermI));

%#####################################################################################################################%
%Calculating the gradient of the likelihood function:
if nargout>1
    gradient(1)=-(sum(I.*(1-IJ).*(f_theta1-(SUMEtheta1./SUME)))-sum(I.*IJ.*(SUMEtheta1./SUME)));
    gradient(2)=-(sum(I.*(1-IJ).*(f_theta2-(SUMEtheta2./SUME)))-sum(I.*IJ.*(SUMEtheta2./SUME))); 
    gradient(3)=-(sum(I.*(1-IJ).*(f_theta3-(SUMEtheta3./SUME)))-sum(I.*IJ.*(SUMEtheta3./SUME))); 
    gradient(4)=-(sum(I.*(1-IJ).*(f_theta4-(SUMEtheta4./SUME)))-sum(I.*IJ.*(SUMEtheta4./SUME))); 
    gradient(5)=-(sum(I.*(1-IJ).*(f_theta5-(SUMEtheta5./SUME)))-sum(I.*IJ.*(SUMEtheta5./SUME))); 
    gradient(6)=-(sum(I.*(1-IJ).*(f_theta6-(SUMEtheta6./SUME)))-sum(I.*IJ.*(SUMEtheta6./SUME))); 
    gradient(7)=-(sum(I.*(1-IJ).*(C-(SUMEtheta7+SUMEInactiveCost.*SUMEInactiveTrips)./SUME))+sum(I.*IJ.*(SUMEInactiveCost-(SUMEtheta7+SUMEInactiveCost.*SUMEInactiveTrips)./SUME))); 
    gradient(8)=-(sum(I.*IJ.*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*(SUMEInactiveTrips./SUME))); 
    gradient(9)=-(sum(I.*IJ.*Y(:,1).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,1).*(SUMEInactiveTrips./SUME))); 
    gradient(10)=-(sum(I.*IJ.*Y(:,2).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,2).*(SUMEInactiveTrips./SUME))); 
    gradient(11)=-(sum(I.*IJ.*Y(:,3).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,3).*(SUMEInactiveTrips./SUME))); 
    gradient(12)=-(sum(I.*IJ.*Y(:,4).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,4).*(SUMEInactiveTrips./SUME))); 
    gradient(13)=-(sum(I.*IJ.*Y(:,5).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,5).*(SUMEInactiveTrips./SUME))); 
    gradient(14)=-(sum(I.*IJ.*Y(:,6).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,6).*(SUMEInactiveTrips./SUME))); 
    gradient(15)=-(sum(I.*IJ.*Y(:,7).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,7).*(SUMEInactiveTrips./SUME))); 
    gradient(16)=-(sum(I.*(1-IJ).*(f_theta16-(SUMEtheta16./SUME)))-sum(I.*IJ.*(SUMEtheta16./SUME))); 
end
